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In this letter we present a real space density functional theory (DFT) localized basis set semi- 
empirical pseudopotential (SEP) approach. The method is applied to iron and magnesium oxide, 
where bulk SEP and local spin density approximation (LSD A) band structure calculations are shown 
to agree within approximately 0.1 eV. Subsequently we investigate the qualitative transferability of 
bulk derived SEPs to Fe/MgO/Fe tunnel junctions. We find that the SEP method is particularly 
well suited to address the tight binding transferability problem because the transferability error at 
the interface can be characterized not only in orbital space (via the interface local density of states) 
but also in real space (via the system potential). To achieve a quantitative parameterization, we 
introduce the notion of ghost semi-empirical pseudopotentials extracted from the first-principles cal- 
culated Fe/MgO bonding interface. Such interface corrections are shown to be particularly necessary 
for barrier widths in the range of 1 nm, where interface states on opposite sides of the barrier cou- 
ple effectively and play a important role in the transmission characteristics. In general the results 
underscore the need for separate tight binding interface and bulk parameter sets when modeling 
conduction through thin heterojunctions on the nanoscale. 



I. INTRODUCTION 

Recently, first-principles theoretical predictions of 
crystalline Fe/MgO/Fe tunneling magnetoresistance 
(TMR) on the order of several hundred percent or mord^ 
were confirmed in series of notable experiments P^l With 
optimization efforts continuing, this dramatic TMR en- 
hancement has placed magnetic tunnel junction (MTJ) 
devices in a unique position to revolutionize memory, 
magnetic sensor, and computing technologies 

The large tunneling magnetoresistance of crystalline 
Fe/MgO/Fe junctions can be understood in terms of 
the symmetry of the MgO crystal, which allows states 
with Ai symmetry to transmit efficiently through the 
band gap of MgO while states of A 2 /5 symmetry decay 
rapidlyfS Near the Fermi energy, the Fe majority and 
minority states are primarily of Ai and A5 symmetry 
respectively. Therefore the MgO barrier acts as a spin 
filter, resulting in half-metallic like conduction between 
Ai states on opposite sides of the barrier. Studies have 
shown that only a single crystalline Fe layer adjacent 
to MgO is sufficient to produce most of the TMR ob- 
served in thicker Fe/MgO/Fe devices!^ Recently, it was 
also suggested that spin torque transfer largely occurs 
at the Fe/MgO interface P Therefore it is essential that 
any Fe/MgO/Fe MTJ device transport model correctly 
capture the physical properties of the Fe/MgO interface. 

From a computational perspective, the scalability of 
density functional theory (DFT) in magnetic metals 
presents serious limitations! 10 ^ 11 ! For example, the study 
of spin torque^ and TMR through large scale MTJ cross 
sections inters persed with magnetic impurities and/or 
crystal defects ^ * 14 * 15 ! would be computationally pro- 



hibitive. Scalability is particularly problematic in non- 
collinear magnetic tunnel junction systems, where the 
calculation convergence time can be prodigiou d 10 * 16 ! due 
to the additional spin degree of freedom. Furthermore, 
the tendency of DFT to underestimate semiconductor 
and insulator band gaps limits its capability to quantita- 
tively model device transport characteristics. For exam- 
ple in magnetic tunnel junctions, the co mmonly a pplied 
local spin density approximation (LSDA p * 16 * 17 * 1 ^ signif- 
icantly underestimates the MgO band gap and there- 
fore over estimates the tunneling current and spin-torque 
transfer. In light of these concerns, we are motivated to 
study the applicability of employing the semi-empirical 
pseudopotential method^ in the context of Fe/MgO/Fe 
magnetic tunnel junctions (MTJs). 

The semi-empirical pseudopotential method (herein 
known by the abbreviation SEP) assumes that the 
Hartree and exchange-correlation potential interaction 
between electrons in a crystal lattice can be well approx- 
imated by an angular dependent or spherically symmet- 
ric potential situated at each atomic site. In its sim- 
plest form, where we assume a spherically symmetric 
SEP, the approach is analogous to the atomic sphere 
approximation The SEP approximation was first ap- 
plied to plane wave calculations and benchmarked with 
respect to the bulk properties of Si and CdSeP^ The im- 
plementation was later scaled up and applied to the study 
of quantum dot systems^ possessing a large number of 
atoms. By optimizing the SEP parameter set one is able 
to correct the band gap of the modeled material while 
maintaining DFT wavefunction accuracy!^ The latter 
feature is of utmost importance in Fe/MgO/Fe tunnel 
junctions, where wavefunction symmetry plays a pivotal 
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role in the device transport characteristics. Furthermore, 
this approximation removes the need for a self-consistent 
convergence loop and therefore allows for the study of 
much larger systems. The metho d is also appealing from 
a tight-binding perspective] 21 * 22 ! since it offers the same 
computational advantages and yet is able to rapidly pro- 
duce an accurate parameterization without employing so- 
phisticated optimization algorithms! 2 ^ 

Building upon previous theoretical studies EmUHH we 
examine the applicability of employing semi-empirical 
pseudopotentials (SEPs)^ for the study of Fe/MgO/Fe 
tunnel junctions within a real space localized basis set 
calculations (rather than plane wave calculations^). The 
discussion is divided in two parts. Firstly, the SEP ex- 
traction method is described in detail. Secondly, we eval- 
uate the SEP method with respect to bulk, interface and 
thin barrier parameterizations. 

The method is first benchmarked against bulk Fe and 
bulk MgO LSDA band structure calculations. Subse- 
quently, we show that the bulk derived SEPs are unable 
to quantitatively capture the LSDA derived Fe/MgO/Fe 
interface and thin barrier tunneling characteristicsPI23 
To overcome this shortcoming we therefore introduce a 
separate interface parameterization through the concept 
of ghost semi-empirical potentials localized between the 
Fe and MgO interface atoms. With these interface cor- 
rections, we are then able to quantitatively capture DFT 
tunneling through thin barriers. It is shown that an ac- 
curate interface parameterization is required for barrier 
widths in the range of 1 nm, where interface states on 
opposite sides of the barrier can couple strongly. We 
also evaluate the transferability and importance of MgO 
barrier band gap corrections with respect to the total 
barrier transmission. In general the results underscore 
the need for separate interface and bulk parameteriza- 
tion sets when modeling electron transport through thin 
tunnel junctions. 



II. METHOD 

We briefly outline our simulation method in this sec- 
tion in two parts. Firstly, we outline the chosen local 
atomic orbital DFT method. Secondly, we discuss the 
real space semi-empirical pseudopotential approximation 
applied in this work. The self-consistent non-equilibrium 
green's function (NEGF) DFT transport method applied 
in this work has been discussed extensively in previous 
publications PH 



A. Local Atomic Orbital DFT Method 

The local atomic orbital pseudopotential DFT time in- 
dependent Hamiltonian can be expressed as, 

H = -\v 2 + V ps (r) + V H (r) + V xc [p(r)}, (1) 



where V ps is the pseudopotential term, V H is the Hartree 
term, V xc is the exchange-correlation potential term 
and p is the system charge density. We may ex- 
pand the pseudopotential expression further into local 
and non-local terms following the Klicnman-Bylandcr 
prescription,^ 

V ps (r) = V^(r) + Vl° c (r) (2) 

N 

= y^-(r) + ^ W (|r-r a |). (3) 

a=l 

where a is the atomic index and r Q is a summation 
taken across the pseudopotentials of each atomic posi- 
tion. However, V p " c (r) is usually long ranged (which re- 
duces the sparsity of the Hamiltonian) and therefore also 
computationally problematic. Thus, we screen V p l ° c {vf^ 
by populating the orbitals of the isolated atom and ar- 
rive at a short ranged neutral atom potential, V NA (r), 
for each atomic species. The preferred local atomic or- 
bital Hamiltonian is then written as 

1 N 

H = --^ + V p n s loc (r) + £ V a NA (\r - r a \) 

a=l 

+ 5V H {v) + V xc [p{v)], (4) 

such that the modified Hartree term is given by 
V 2 8V H (r) = -AnSp(r). We define Sp(r) = p(r) - 
TliaP < a 0m ( r ) where p^ om is the neutral atom charge ar- 
rived at by populating the orbitals of an atomic species. 



B. Semi-empirical Pseudopotentials 

1. Extracting SEPs from real space DFT calculations 

In its most basic form, the semi-empirical pseudopo- 
tential approximation^ assumes that all local terms may 
approximated by a spherically symmetric local potential 
around each atom. This is objective is partially accom- 
plished by including V NA (r) but to arrive at a proper 
spherical potential at each atomic site we must also re- 
duce SVh and Vxc such that, 

N 

SV H (r) + V xc [p(r)} « ^d r " r «D- ( 5 ) 

a=l 

where v a is the spherical approximation to the self- 
consistent Hartree and exchange-correlation terms for 
atom a. The full semi-empirical pseudopotential is given 
by V£ EP (r) — v a (r) + V^ A {r). The approach is similar 
in spirit to the atomic sphere approximation applied in 
the muffin-tin orbital method.^ Although not done here, 
angular dependence may be introduced to the SEP term. 
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This leads to a revised Hamiltonian operator, 



H 



SEP 



-V 2 + < oc (r) + 



N 

E^ 5£P ( 

a=l 



.1) (6) 



which does not require a self-consistent loop to solve 
since there is no interdependence between the semi- 
empirical pseudopotentials and the charge density. The 
term "semi-empirical" is applied because these potentials 
are initially derived from first-principles calculations and 
then fitted to experimental data if required - in this work 
to overcome the LSDA band gap underestimation. 
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where V£ denotes a Bessel integral over 8V H (r) + 
V xc [p(r)} and S^ n denotes an overlap integral on the 
right had side of Eq. The coefficients c™ are then 
directly arrived at by matrix inversion. 



2. Extension to bulk systems 



Finally, we would like to extract the SEP for each 
atomic species from a self-consistent DFT calculation. 
Let us assume that the spherical approximation to the 
self-consistent Hartree and exchange-correlation terms, 
v a (r), goes to zero beyond a cutoff radius of r c - which is 
not necessarily equivalent to the cutoff radius of Vj? A (r). 
Note that the zero potential condition outside the cutoff 
radius may be adjusted, for example by adding a positive 
offset to the real space DFT potential. Within the cutoff 
radius we may define a complete orthonormal basi; 
represent v a (r) such that 
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<t>n(r) 



sin(mrr /r c ) 



r <r c 
r > r c 
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v a (r) 
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n=l 



(7) 



(8) 



where the potential is represented by a linear expan- 
sion of the zeroth order spherical Bessel function - the 
eigcnfunctions of an electron with no angular momentum 
trapped in an infinite spherical well of radius r c . Note 
that higher order spherical Bessel functions are not able 
to capture a non-zero system potential at the atomic ori- 
gin. To solve for the coefficients c" we substitute Eq. Q 
into Eq. ^ and construct a linear equation, with N x M 
unknowns, by integrating both sides through with </> m 
centered at atomic species (3, 



|r 

N N 



v \){5V H {v) + V xc [ P {v)])dv 



EE c »/' 



|r-rg|)</> n (|r-r Q |)dr. 



(9) 



a— 1 n— 1 



In Eq. we have forced an equality between the self- 
consistent DFT local potential and the spherical SEP 
approximation to that potential. By further considering 
all Bessel functions in our SEP expansion we obtain linear 



In bulk periodic systems, we need to consider the peri- 
odicity of the system potential when solving for the SEP 
Bessel coefficients. If we take the left hand side of Eq. |9]) 
and integrate through the SEP Bessel functions of a given 
unit cell, the periodic potential (5V H (r) + V xc [p(r)]) 
over which the integral is performed will have contribu- 
tions not only due to the SEPs of the unit cell which 
we have selected but also due to the SEPs of neighbor- 
ing unit cells. We can address this issue by adopting 
a supercell tight binding description of bulk periodicity, 
where beyond twice the maximum SEP cutoff radius the 
interaction between a unit cell and its bulk neighbors is 
assumed to go to zero. In this manner, the SEP Bessel 



integrals on the left hand side of Eq. (10 1 are performed 
only for the central unit cell in our supercell. However, 
the SEP coefficients must be the same for all unit cells. 
Therefore, the SEP matrix overlap matrix on the right 
hand side of Eq. (10 1 is expanded into a summation of 



the SEP overlap matrices between the central unit cell 
and all neighboring unit cells within the supercell. The 
revised unit cell SEP integral equation is then written as, 



•Ml* -t p \){SV b {t) + V xc [p(r)})dv 



N N . 

= E E c ™ E / ^»(l r " r /3l)<M|r - r a - R\)dr. 

1 1 T» J 



(ii) 



such that R = uiR.i + Ti2R2 + ^3R-3; where R-1,2,3 are the 
translation vectors of the unit cell and ri\ 2,3 are integers. 



3. Extension to collinear spin polarized systems 

Thus far we have only outlined the SEP extraction 
procedure for spin independent calculations. When mod- 
eling collinear spin polarized systems, separate SEPs 
are extracted for the majority and minority spin elec- 
trons. For majority spin up electrons we simply set 
SV H + V xc T[p(r)} a £a=i ^(|r - r Q |) in Eq. §, and 
solve for the d^} zeroth order Bessel coefficients following 
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the above prescription. In the same manner one is able 
to extract the minority spin down SEP coefficients. 



4- Band structure fitting methodology 

To fit the SEP band structure of a given bulk mate- 
rial, the Jacobian matrix [J] of target band structure 
points pj is computed with respect to the material SEP 
coefficients c, such that [J]ij = dpj/dci. For example, 
to correct the band gap of insulators or semiconductors, 
the valance and conduction band energies at symmetry 
fc-points can be taken as targets to be raised or lowered. 
After computing the Jacobian, a new set of coefficients 
is computed via 

c new =c + [J]-l p (12) 

where [J] -1 is the pseudo-inverse of the Jacobian if the 
matrix is not square and p is the vector between the ex- 
isting band structure target values (derived from c) and 
the desired band structure target values. The process is 
iterated by setting c = c new and recalculating the Ja- 
cobian, until the vector p approaches a small tolerance 
value (say 0.1 eV per target value). 
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III. SEP ELECTRONIC STRUCTURE FOR 
IRON AND MAGNESIUM OXIDE 

A. Bulk Fe and MgO 

To insure quantitative transport calculations, in agree- 
ment with existing DFT methods, the semi-empirical 
pseudopotential approximation must be benchmarked 
against self-consistent results. Therefore, we begin by 
examining the accuracy of the semi-empirical pseudopo- 
tential method detailed in Sec. [Il]as applied to bulk iron 
and bulk magnesium oxide LSDA calculations. 

The band structure of bulk Fe in the (001) direc- 
tion is presented in Fig. [l^, where (001) is the direc- 
tion of electron transport through Fe/MgO/Fe tunneling 
barriers.E^EIafil The lattice constant of Fe is set at 2.87 
AGHland a } on g ran g e double-^ polarized basis setP^is em- 
ployed in all calculations. The SEP cut off radius is set to 
5 Bohr. The LSDA calculated band structure (presented 
as a dashed blue line) in Fig. [T^ can be seen to agree quite 
well with the SEP calculated band structure (presented 
as a green line) . The mean margin of error between the 
two band structure calculations is approximately 0.1 eV. 

The band structure of strained bulk MgO is presented 
Fig. The MgO lattice consta nt is set at 4.21 A in 
the (001) transport direction) 31 ' 32 1 However the (100) and 
(010) directions are strained by 3.8 %, to 4.05 A, in order 
to lattice match bulk Fe (see the two probe Fe/MgO/Fe 
IIIB[ ). A double-^ polarized basis 



calculations in Sec. 



seP^I is employed in all calculations. The Mg atoms are 
assigned a basis set a cut off radius of 8 Bohr and the 



FIG. 1: Bulk band structure of Fe and MgO. Subfigure a) pro- 
vides the Fe (100) crystal bulk band structure with the Fermi 
energy is situated at eV. Subfigure b) provides the strained 
MgO bulk band structure with the Fermi energy positioned to 
match that of MgO sandwiched between two Fe(100) slabs. 
The LSDA calculated band structure is shown as a dashed 
blue line. The LSDA SEP fit is shown in as a solid green line. 
The modified LSDA SEP result fitted to the bulk MgO band 
gap of 7.7 eV^is shown as a dot-dashed gold line. 



O atoms a cutoff radius of 4.5 Bohr. The LSDA calcu- 
lated band structure (presented as a dashed blue line) in 
Fig. [TJd can be seen to agree quite well with the SEP cal- 
culated band structure (presented as a solid green line). 
The solid green band structure in Fig. [TJd imposes SEP 
cutoff radii of 5 Bohr and 4.5 Bohr to Mg and O re- 
spectively. The margin of error between the LSDA and 
SEP band structures calculations is approximately 0.1 
eV. We have found the same level of SEP fit accuracy 
can be achieved with the unstrained MgO lattice. 

It is important to note that SEPs with longer cutoff 
radii (beyond 5 Bohr as demonstrated here) are tenable 
but often end up sampling not only the potential of the lo- 
cal atom which they are situated on but also the potential 
of neighboring atoms. Such SEPs are therefore not even 
qualitatively transferable to material heteroj unctions (for 
example Fe/MgO/Fe as studied in this work). 

To achieve the above fit accuracy, with both bulk Fe 
and MgO, we applied a LSDA real space grid resolution 
of 4 points per Bohr (64 points per Bohr 3 ). To fit the 
Fe LSDA band structure, 20 Bessel functions for both 
the up-spin and down-spin SEPs were required, although 
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with shorter range Fe basis sets we have found that as few 
as 10 Bessel functions are suitable. To fit the MgO band 
structure 10 Bessel functions per SEP were required. Re- 
ducing real space the grid resolution reduces the accuracy 
of the integrals in Eq. ^ and can result in a poor match- 
ing between the LSDA and SEP calculated band struc- 
tures. Likewise, an insufficient number of Bessel func- 
tions in Eq. Q will result in a poorly constructed SEP. 
There is a fine balance between the grid resolution and 
the number of Bessel functions, as too much of either can 
raise both the calculation computation time and memory 
consumption. 

Radial real space plots of the strained bulk MgO and 
bulk Fc SEPs are presented in Fig. |2j The O SEP is 
much sharper than the Mg SEP (see solid green lines in 
Figs. [2^, and read off the right axis), and both pos- 
sess considerable corrections when the MgO band gap 
is expanded (see dot-dashed gold lines in Figs. [2^i and 
[2]d read of the left axis). The SEP corrected MgO band 
structure, fitted to the bulk MgO band gap of 7.7 eVP 
is plotted as dot dashed gold line in Fig. [I]}. We have in- 
vestigated shorter ranged band gap corrections but have 
found that they are not able to open the band gap with- 
out significantly distorting the band structure. The Fe 
SEPs are displayed in Fig. [2ji. The up spin Fe SEP (see 
solid green line in Fig. |2j; read off the right axis) and the 
down spin Fe SEP (see double-dot-dashed black line in 
Fig. [2J: read off the left axis) differ primarily only with 
respect to on-site exchange corrections localized at the 
Fe atomic core. 
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FIG. 2: Spherical SEPs for each of the elements. Subfigures 
a) and b) show the LSDA Mg and O SEPs for MgO in solid 
green and the MgO band gap fit correction to the LSDA SEPs 
in dot-dashed gold. The MgO band gap fit corrections should 
be read off the left axis and the LSDA SEPs should be read off 
the right axis. Subfigure c) shows the LSDA Fe up spin SEP 
in solid green (read off the right axis) and difference between 
the Fe down spin and Fe up spin LSDA SEPs {V§^ p - V§^ p ) 
as a double-dot-dashed black line (read off the left axis). The 
vacuum level is set at eV. 



B. Bulk SEP Transferability to Fe/MgO/Fe Tunnel 
Junctions 

Given the accurate bulk band structure results pre- 
sented in the previous section, we now proceed to ex- 
amine the transferability of bulk derived semi-empirical 
pseudopotentials to magnetic tunnel junctions. The 
bulk SEP zero bias two probe Fe/MgO/Fe TMR ra- 
tio, projected density of states (PDOS), and transmis- 
sion characteristics are shown to qualitatively match 
first-principles self-consistent NEGF-LSDA results. 1121251 
To obtain a quantitative tunneling barrier parameter- 
ization we explicitly identify the real space Hartree 
and exchange-correlation potential error of bulk SEPs 
at the Fe/MgO interface, and introduce the notion of 
ghost SEPs to fit and thereby remove the transferabil- 
ity error. In this regard the SEP tight binding ap- 
proach is shown to be advantageous as it allows both 
orbital space and real space characterization of het- 
croj unction interface errors introduced by bulk param- 
eterizati ons. N on-pseudopotential based tight binding 
methods] 21 * 23 * 34 ^ where the atomic orbitals overlap inte- 
grals are used as fitting parameters, do not allow such 
a systematic characterization of interface transferability 
errors. 



1. Fe/MgO/Fe MTJ Geometry 



The Fe/MgO/Fe tunnel junction under investigation 
consists of 5 MgO layers P The full NEGF-LSDA device 
region is shown in Fig. [3] where the semi- infinite leads are 
accounted for by self-energy terms in the device Green's 
function.^ We have set the Fe lattice constant in both 
leads to 2.87 A - the MgO transverse lattice constant is 
also set at 4.05 A. However, the MgO layers are sepa- 
rated by 2.1 A in the transport di rectio n, matching the 
bulk MgO lattice constant of 4.2 A. 13JJ23J The Fe-0 bond- 
ing distance at the Fe/MgO interface is set at 2.169 A.- 
The unit cell geometry shown in Fig. [3] is periodically 
repeated infinitely in the transverse (x, j/)-plane (which 
lies perpendicular to the tunneling transport z-direction) . 
We have chosen to examine the five layer MgO device ge- 
ometry, rather than wider or thinner barriers, because at 
this thickness the bulk MgO band gap reappears in the 
middle of the barrier. This allows a proper evaluation of 
both the interface and bulk properties of the MgO tun- 
neling barrier as approximated by the SEP method. 
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FIG. 3: Fe/MgO/Fe 5 layer device geometry. The atomic 
color index is as follows: iron atoms are colored gold, the 
magnesium atoms are colored green, and oxygen atoms are 
colored red. The system is mirror symmetric along the z-axis 
about layer 11 (the middle of the barrier). 



2. Fe/MgO/Fe MTJ Potential Profile Study 

Given that our SEP approach relies upon the funda- 
mental assumption that the potential of a system can 
be approximated by a summation of local potentials, 
we begin by comparing the system potential results of 
two probe SEP and self-consistent NEGF-LSDA calcu- 
lations. This SEP tight binding method allows not only 
orbital space evaluation, in the form of projected density 
of states (PDOS) plots, but more importantly real space 
evaluation of the tight binding Hamiltonian. 

In Figs. [4^ and |4Jd potential cuts, through the Mg and 
O interface atoms respectively, of the MgO two probe ge- 
ometry (see Fig. [3]) are plotted in the electron transport 
z-direction. The total down spin potential of our 5 layer 
Fe/MgO/Fe device is given as a dotted black line (to be 
read off the right axis) and the corresponding two probe 
bulk SEP transferability error is plotted in solid green 
(to be read off the left axis) . The up spin and down spin 
bulk SEP transferability errors are very similar, therefore 
in the interest of a concise discussion we include only the 
down spin results. Lastly, the MgO SEP potentials dis- 
cussed in this section do not include a band gap cor rec- 
tion (see Figs. [I] and [2]), this issue left to Sec. IIIB4 



Away from the interface the SEP potential error, 
shown in green in Fig. [4] and read off the left axis, is 
largely flat apart from small oscillations on the Fe and 
Mg atoms and peaks localized on the O atoms. The small 
oscillations away from the interface, can be attributed to 
the spherical approximation where we have neglected an- 
gular variations in the crystal potential about an atom. 
The sharp errors localized on each oxygen atom are due 
the small number of Bessel functions (ten per atom) em- 
ployed in the bulk MgO fit, which are not able to com- 
pletely capture the rapid drop in the system potential at 
the oxygen atomic core. However, due to their sharp na- 
ture these peaks contribute negligibly to the integrated 
Hamiltonian oxygen onsitc energies and therefore can be 
ignored (see Sec. IIIB3 for further details). Immediately 
away from the interface, the bulk and two probe system 
potentials agree remarkably well. However, at the inter- 
face the SEP potential error is substantial. It is impor- 
tant to note that bulk MgO SEP and LSDA two-probe 
MgO Fermi energies have been aligned via a constant 
bulk potential shift (see Fig. [T}d) . 

The interface bulk SEP error is localized at the FeMg 
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FIG. 4: The total NEGF-LSDA spin down potential (V To t = 
^ r Tot DA ) i s shown as a dotted black line, with its axis given 
on the right hand side of the figure (the vacuum level is set 
at eV). The bulk LSDA SEP interface error is shown in 
green (AV To t = VSM P ~ VTot DA ) and the interface spherical 
SEP correction to the error is shown in red (both potentials 
are read off the left axis). Subfigure a) displays the system 
potential as a linear cut in the z-direction through the Mg 
atom at the Fe/MgO interface. Subfigure b) displays the sys- 
tem potential as a linear cut in the ^-direction through the 
FeO bond at the Fe/MgO interface. An atomistic cartoon is 
shown to scale above each potential plot, where a dip in the 
total potential corresponds to an atomic nuclear position - Fe 
atoms are gold, O atoms are red, and Mg atoms are green. 
Up spin results are nearly identical. The SEP MgO band 
gap corrections (see Figs, [I] and [5| are not included in this 
comparison. 



junction, see the solid green line Fig. [3Ji, and reaches a 
maximum of approximately 4 eV (read off the left axis of 
Fig. |4ji) . Yet, the integ rated local atomic orbital matrix 
Hamiltonian errors^ESl occur on both the Mg interface 
atoms and the oxygen bonded Fe interface atoms. On 
the other hand, the FeO potential cut (in Fig. |4jo) dis- 
plays relatively little error - although this error is slightly 
larger for the up spin system potential. By including 
the bulk Fe and bulk MgO Fermi level energy offset in 
our bulk MgO SEP fit, we have largely compensated for 
the Hartree potential created by charge redistribution at 
the FeO interface^ (a classical analogue to this would 
be the built in potential profile of a semiconductor p-n 
junction). This offset minimizes the FeO bonding poten- 
tial error, which can be largely attributed to neglected 
changes in the exchange-correlation potential. However, 
at the FeMg interface the bonding environment changes 
even more drastically, each Mg interface atom loses one 
nearest neighbor, and the charge redistribution cannot be 
approximated by the Hartree potential required to align 
the heterojunction Fermi energies. By removing a near- 
est neighbor at the FeMg interface we violate the spheri- 
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cal symmetry that our SEP fit assumes for the chemical 
environment and therefore a fundamentally asymmetric 
solution to the SEP approximation is required to over- 
come the transferability error. 

To overcome the bulk transferability error we introduce 
the concept of ghost SEPs, that is SEPs which are not 
localized at an atomic core but instead situated within 
the bonding region of the heterojunction interface. Such 
ghost SEPs are fitted to cancel the SEP interface trans- 
ferability error, that is the potential difference between 
the LSDA two probe calculation and the bulk SEP two 
probe calculation. In this manner ghost SEPs allow sep- 
arate DFT bulk and heterojunction interface paramctcr- 
izations, which can be applied independently (for exam- 
ple) to study the device transport properties of various 
barrier widths, spin torque,^ and the role of electron 
and spin defect/impurity scattering within the barrier.^ 
The interface bulk SEP transferability error is analyzed 
in further detail in Sec. |III B 3| and Sec. |III B 4| with re- 
spect to the two probe interface PDOS and transmission 
properties. 

The down spin ghost SEP parameterization applied 
throughout this work is shown as a solid red line (read 
off the left axis) in Figs, [iji and |4j). We have employed 
two spherically symmetric ghost SEPs per spin (four to- 
tal) localized along the FeO and FeMg line cuts as shown 
in Fig. [4j although higher order angular momentum SEPs 
may also be applied to capture asymmetry at the inter- 
face. The SEPs plotted in Fig. [3] posses a cutoff radius 
of 5 Bohr and are composed of 10 Bessel functions. The 
Bessel coefficients of the ghost SEPs are arrived at by 
replacing the terms SV H (r) + V xc (r) in Eq. 1 1 with the 
two probe potential difference V^ t DA (r) 



Bulk Fit 



Ghost Fit 



3. Fe/MgO/Fe MTJ PDOS Transferability Study 

To conceptualize the importance of a proper interface 
parameterization, let us begin by comparing the PDOS 
results of two probe SEP and self-consistent NEGF- 
LSDA calculations. 

The parallel orientation PDOS results at layers 6, 8, 9 
and 11 in the Fe/MgO/Fe device are displayed in Fig. [EJ 
- see the geometry diagram in Fig. [3] for details on the 
layer numbering. The two probe NEGF-LSDA PDOS is 
shown in dashed blue, the bulk SEP PDOS in green, and 
the ghost SEP PDOS in solid red. It is important to note 
the mirror symmetry of the 5 layer Fe/MgO/Fe system, 
where under zero bias conditions the PDOS at layers 6, 
8 and 9 is equivalent to the PDOS at layers 16, 14 and 
13 respectively. 

As we transition from deep within the Fe leads towards 
the MgO tunnel junction, the bulk SEP, NEGF-LSDA, 
and ghost SEP PDOS calculations agree well. This agree- 
ment holds up until the third Fe layer as measured from 
the Fe/MgO interface - see the layer 6 PDOS results in 
Fig. [5j The disagreement reaches a maximum directly at 
the Fe/MgO interface (see the bulk SEP and ghost SEP 
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FIG. 5: Fe/MgO/Fe SEP parallel orientation projected den- 
sity of states (PDOS) before approaching the interface (layer 
5), at the interface (layers 8 and 9), and in the middle of the 
barrier (layer 11) - see Fig. [3] for details on the layer num- 
bering. The Fermi energy is located at eV, the imaginary 
green's function broadening is set at 25 meV, and the fc-point 
sampling set at 8x8 in the transverse Brillouin zone. The up 
spin PDOS is shown on the positive axis and the down spin 
PDOS is shown on the negative axis. The reference two probe 
NEGF-LSDA PDOS is shown as a dashed dark blue line in 
each figure. The bulk SEP fit (displayed in solid green) ap- 
plied to the two probe geometry is shown in the first column. 
When ghost SEPs (displayed in solid red) are introduced at 
the interface to correct the bulk Fe and MgO SEP trans- 
ferability errors, an accurate fit is obtained as shown in the 
second column. 



fits respectively at layers 8 and 9 in Fig. [EJ where the 
bulk SEP PDOS begins to diverge from the LSDA PDOS. 
Further within the MgO barrier (layer 11 in Fig. IE) the 
bulk SEP and ghost SEP parameterizations both cap- 
ture the LSDA calculated MgO band gap as it begins to 
reappear. The agreement between the SEP and LSDA 
results at layer 11 Fig. [5] clearly shows that that sharp 
SEP fit errors located on the O atoms in Fig.|4]do not in- 
fluence the barrier electronic structure (see the discussion 
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in Sec. IIIB2) 



Returning to the interface, we see that the layer 8 
NEGF-LSDA result (in dashed dark blue in Fig. |>} dis- 
plays the characteristic Fe/MgO interface PDOS includ- 
ing the minority PDOS Fermi energy resonant peak! 1 ! 36 ! 
Yet if we turn our attention to the bulk SEP PDOS layer 
8 result (shown in green in Fig. |5| we see a noticeable 
disagreement, namely the characteristic twin peak (A 5 
and Ai as labeled in Fig. [5]) Fe minority and majority 
PDOS resonances are markedly distorted and in the case 
of the majority interface states there is a further 0.75 eV 
upwards shift. This upwards shift in the interface bulk 
SEP majority interface states is due entirely to the pos- 
itive nature of the interface error as shown by the green 
potential plotted in Fig. [3Jl 

The Ai interface state decays slowly into the MgO 
barrier and the A5 interface state decays rapidly into 
the MgO barrier (see Sec. IIIB4). Furthermore, the 2 



eV exchange splitting between the majority and minority 
carriers results in half-metallic like conduction between 
the slowly decaying Ai interface states, which dominate 
the TMR and spin torque characteristics of Fe/MgO/Fc 
junctions!^ By first distorting the Ai minority/majority 
interface states and then shifting the majority Ai inter- 
face state by 0.75 eV, the bulk SEP approximation in- 
troduces considerable error into the half-metallic proper- 
ties of the the Fe/MgO/Fe tunneling as we show in fur- 
ther detail in the next section. However, the ghost SEPs 
clearly (as shown in red in the second column of Fig. [5| 
are able to almost entirely compensate the for bulk SEP 
PDOS interface transferability errors. This quantitative 
result is achieved with only first order (spherically sym- 
metric) ghost SEPs, where angular dependent interface 
ghost SEPs might be necessary for more complex hetero- 
junction interfaces. 



4. Fe/MgO/Fe MTJ Transmission Study 

Thus far we have performed a detailed analysis of 
the interface potential and PDOS errors (see Figs. [4] 
and [5]) which result when bulk SEPs are transfered to 
Fe/MgO/Fe tunnel junctions. Yet, for the purposes of 
electron device modeling we are most interested in the 
transport implications of such interface errors. In this 
regard, previous studies have shown resonant i nterface 
states can signi fica ntly influence the transmission 31 ! 36 ! 37 ! 
and spin torque^^S! properties of MTJ barriers. 

The zero-bias total transmission of our Fe/MgO/Fe 
tunneling device geometry (see Fig. [3| is presented in 
Fig. [6] The LSDA transmission is shown in dashed 
dark blue, the bulk SEP transmission in green, the ghost 
SEP transmission in solid red, and the MgO band gap 
corrected transmission is shown in dotted black (recall 
the dot-dashed gold correction potentials displayed in 
Figs. [2^i and[2|)). The zero bias parallel orientation up 
spin transmission is displayed in Fig. [6^,, the parallel ori- 
entation down spin transmission is displayed in Fig. [6b, 



and the antiparallel transmission is displayed in Fig. [6J;. 

An initial inspection of Figs. [6^, through [6]:, both above 
and below the Fermi energy where the Fe/MgO interface 
states are more prominent, reveals sizable transmission 
corrections between the bulk SEP result (in solid green) 
and the ghost SEP result (in solid red) - where the latter 
matches the LDSA transmission (in dashed blue) quan- 
titatively. The Fe/MgO/Fe transmission characteristics 
are determined by the rapid decay of the A5 interface 
state and slow decay of the Ai interface state within 
MgO, resulting in half-metallic like tunneling between 
majority and minority Ai Fe/MgO interface states on 
opposite sides of the barrier (see layer 8 in Fig. [jjl. In 
Fig.[6]the half-metallic like conduction is immediately ev- 
ident, where at -1 eV in Fig. |6|i and at +1 eV in Fig. [BJa 
we see a large rise in transmission corresponding to the 
onset of the majority and minority spin Ai states respec- 
tively (also shown on layer 8 in Fig. [5]) . 

The Fermi energy (set at eV in Figs. [5] and [6]) low bias 
TMR for the LSDA, bulk SEP and ghost SEP methods 
is calculated to be 99.85, 59.39, 97.47 respectively. The 
bulk SEP approximation underestimates the Fe/MgO/Fe 
low bias TMR by 30%, and the ghost SEP interface cor- 
rections are able to compensate quite accurately. This 
low bias error is due entirely to the Ai majority interface 
state broadening/shift and as shown on layer 8 in Fig. [5] 
and in Fig. [6^ where the bulk SEP (solid green) transmis- 
sion is shifted by 0.5 eV. The minority spin A 5 interface 
state is not significantly altered by the bulk SEP approx- 
imation (see layer 8 in Fig. [5]) resulting in little change 
in the half-metallic like transmission at the Fermi energy 
of the Ax interface states in the anti-parallel orientation 
(see Fig. |6j;). The same holds for the parallel minority 
transmission although it's contribution to the low bias 
parallel current is negligible (see Fig. (6j}). 

Though the minority spin A5 interface state transmis- 
sion is not significantly altered by the bulk SEP (solid 
green) transferability error, the minority spin Ai inter- 
face state zero bias transmission is however drastically 
underestimated between eV and 1 eV (see Fig. [HJd) . 
Similarly, we can see a notable underestimation of the 
majority spin bulk SEP (solid green) Ax zero bias trans- 
mission between -1 eV and eV as shown in Fig. [6ji. 
Given the half-metallic like spin filtering property of the 
MgO barrier, in which tunneling between Ai states domi- 
nates, these errors have a significant impact on the biased 
tunneling current. Without the ghost SEP interface cor- 
rections, the parallel tunneling current is underestimated 
by a factor of 2 and the antiparallel current by up to an 
order of magnitude within the bias window of 1 V as 
shown in Fig. [6ji (compare the ghost SEP solid red and 
bulk SEP solid green results). Furthermore, in Fig. [6)5 the 
bulk SEP (solid green) Fe/MgO/Fe TMR displays sizable 
deviations from the interface corrected ghost SEP (solid 
red) TMR, lacking the characteristic smooth decay under 
biasP 

Similar, though less sizable, interface errors occur 
when we introduce bulk SEP MgO band gap correc- 
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FIG. 6: Subfigures a), b) and c) display the parallel and anti-parallel zero bias total transmission with respect to energy through 
the 5 layer Fe/MgO/Fe MTJ geometry shown in Fig. [3] The self-consistent LSDA total transmission is shown in dashed dark 
blue, the bulk SEP result is shown in solid green, the ghost SEP result is shown in solid red, and the band gap corrected 
ghost SEP result is shown in dotted black. The biased device parallel (upward pointing triangles) and anti-parallel (downward 
pointing triangles) currents for the ghost SEP (red solid line) and bulk SEP (green solid line) approximations are shown in 
subfigure d) in units of Amperes per unit cell (2.87 A by 2.87 A). The voltage profile is assumed to drop linearly across the 
MgO barrier for the calculated IV points. The ghost SEP (red solid line) and bulk SEP (green solid line) TMR ratios under 
bias are shown in subfigure e). The Fermi energy is set at eV, the imaginary green's function broadening is set at 25 meV, 
and the fc-point sampling set at 100x100 in the transverse Brillouin zone. We define TMR — (Ip — Iap)/Iap, where Ip is the 
parallel current and Iap is the antiparallel current. 



tions to our ghost SEP Fe/MgO/Fe Hamiltonian (sec 
dotted black transmission plots in Figs. [6^1 through^). 
The bulk SEP MgO band gap corrections are plotted in 
Figs. |2^l and Figs. (2Jd (see gold dot-dashed lines read off 
the left axis). Reaching up to 5 Bohr, these band gap 
corrections extend into the Fe/MgO interface and suffer 
the same transferability problem as the uncorrected bulk 
SEPs. However, from Figs. [6^ through [6ji it is evident 
that the primary role of the band gap correction is to 
lower the tunneling current by an order of magnitude (as 
expected). Likewise, the band gap correct TMR ratio at 
43.94 compared to the LSDA TMR ratio of 99.85, can 
be attributed to the reintroduced interface state errors 
rather than a fundamental alteration in the Fe spin fil- 
tering properties of MgO P However, it may be necessary 
in future studies to simultaneously address the nature of 
Fe/MgO interface states and MgO exchange-correlation 
corrections beyond the LSDA approximation.^ 

IV. SUMMARY 

We have detailed a straight forward method for ex- 
tracting semi-empirical pseudopotcntials from real space 
DFT calculations. The method has been shown to pro- 
duce accurate bulk derived spherical SEPs, matching self- 
consistent LSDA band structure results to within 0.1 eV. 
Subsequently, we examined the transferability of bulk de- 
rived MgO and Fe SEPs to Fe/MgO/Fe tunnel junctions. 
It was shown that LSDA calculated Fe/MgO interface 
states are not adequately described by bulk SEPs. As 
a result bulk SEPs can significantly underestimate or 
overestimate of the spin dependent transmission through 



thin Fe/MgO/Fe tunnel junctions. However, the SEP 
tight binding method allows characterization of both the 
system layer by layer PDOS and real space potential. 
Primarily due to the local chemical environment change 
experienced by interface Mg atoms, where the number 
of nearest neighbors is reduced from six to five, an in- 
terface fit is required to overcome bulk transferability 
errors. Therefore, we put forward the notion of ghost 
SEPs, not localized at an atomic site but within the 
interface, to parameterize DFT calculated heterojunc- 
tions. The ghost SEPs interface parameterization was 
shown to completely recover the LSDA interface PDOS 
and Fe/MgO/Fe transmission characteristics. In general, 
the results emphasize the need for separate bulk and in- 
terface parameterizations when applying tight binding 
methods to study transport through nanoscale hetero- 
junctions where interface states can couple effectively. 
Lastly, we note that all the parameters and the device 
geometry discussed herein are provided for download as 
supplementary materialP*l 
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